/*
 * PCA.cpp
 *
 *  Created on: 10.11.2009
 *  Author: 	Martin Do
 *	Copyright: 	Martin Do, Chair Prof. Dillmann (IAIM),
 *            	Institute for Anthropomatics (IFA),
 *		        Karlsruhe Institute of Technology (KIT). All rights reserved.
 *
 *	Changes:
 *
 *	Description:
 *
 *	Mit Eigenvektoren oder mit EM
 */



#include "PCA.h"
#include <math.h>
#include <assert.h>

CPCA::CPCA()
{

}

CPCA::~CPCA()
{
}

void CPCA::computePCA_Cov(std::vector<gsl::vector> dataSet, int nReducedDimensions)
{
	
	matrix datamatrix;
	
	datamatrix = CMatrixOperation::getMatrixFromData(dataSet);
	
	
/************************
Implement PCA
Calculate mProjectionMatrix (member)
   ***********************/
   
   matrix eigenvectors;
   
   mRowMeans = CMatrixOperation::getMeanVector(datamatrix);
   mCovMatrix = CMatrixOperation::getCovarianceMatrix(datamatrix, mRowMeans);
   
   eigenvectors = CMatrixOperation::getEigenValues(mCovMatrix, true);
   mProjectionMatrix = matrix (dataSet[0].size(), nReducedDimensions);
   
   //mProjectionMatrix = gsl_matrix_alloc (dataSet[0].size(), nReducedDimensions);
   
   for (unsigned int i = 0; i < dataSet[0].size(); i++) {
       for (int j = 0; j < nReducedDimensions; j++) {
		   mProjectionMatrix(i,j) = eigenvectors(i+1,j);
	   }
   }
	

}

std::vector<gsl::vector> CPCA::reduceDim(std::vector<gsl::vector> dataSet)
{
	matrix dataMatrix = CMatrixOperation::getMatrixFromData(dataSet);
	
    for (unsigned int i = 0; i < dataMatrix.get_cols(); i++)
            for (unsigned int j = 0; j < dataMatrix.get_rows(); j++)
				dataMatrix(j,i) = dataMatrix(j,i)-mRowMeans(j);
	
	matrix reducedDataMatrix = dataMatrix.transpose() * mProjectionMatrix;
	
	
  /************************
Implement projection to feature space using mProjectionMatrix
   ***********************/
   std::vector<gsl::vector> reducedDataSet = CMatrixOperation::getDataFromMatrix(reducedDataMatrix);	
   return reducedDataSet;
}

std::vector<gsl::vector> CPCA::restoreDim(std::vector<gsl::vector> dataSet)
{
	matrix dataMatrix = CMatrixOperation::getMatrixFromData(dataSet);
	matrix restoredDataMatrix = mProjectionMatrix * dataMatrix.transpose();
    for (unsigned int i = 0; i < restoredDataMatrix.get_cols(); i++)
            for (unsigned int j = 0; j < restoredDataMatrix.get_rows(); j++)
				restoredDataMatrix(j,i) = restoredDataMatrix(j,i)+mRowMeans(j);
	std::vector<gsl::vector> restoredDataSet = CMatrixOperation::getDataFromMatrix(restoredDataMatrix);
	return restoredDataSet;
}
